Viscoelastic Phase Separation in Shear Flow 
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We numerically investigate viscoelastic phase separation in polymer solutions under shear using a 
time-dependent Ginzburg-Landau model. The gross variables in our model are the polymer volume 
fraction and a conformation tensor. The latter represents chain deformations and relaxes slowly 
on the rheological time giving rise to a large viscoelastic stress. The polymer and the solvent obey 
two-fluid dynamics in which the viscoelastic stress acts asymmetrically on the polymer and, as a 
result, the stress and the diffusion are dynamically coupled. Below the coexistence curve, interfaces 
appear with increasing the quench depth and the solvent regions act as a lubricant. In these cases 
the composition heterogeneity causes more enhanced viscoelastic heterogeneity and the macroscopic 
stress is decreased at fixed applied shear rate. We find steady two-phase states composed of the 
polymer-rich and solvent-rich regions, where the characteristic domain size is inversely proportional 
to the average shear stress for various shear rates. The deviatoric stress components exhibit large 
temporal fluctuations. The normal stress difference can take negative values transiently at weak 
shear. 

PACS numbers: PACS numbers: 64.75. +g, 83.60.Df, 83.60.Rs, 83.80.Rs 



I. INTRODUCTION 

In phase-separating polymer systems, domain mor- 
phologies are influenced by a number of factors, including 
the molecular weights, the composition, closeness to the 
critical point, and viscoelasticity 0,13 • Particularly when 
the two components have distinctly different viscoelastic 
properties as in semidilute polymer solutions and poly- 
mer blends of long and short chains, unique interplays 
emerge between viscoelasticity and thermodynamic in- 
stability. In such asymmetric systems, salient effects ob- 
served in experiments are as follows. First, in early stage 
spinodal decomposition, the kinetic coefficient L(q) de- 
pends on the wave number q as L(q)/L(0) ~ (Cve?) -2 for 
q larger than the inverse of a viscoelastic length £ ve 
where £ V c can be much longer than the gyration radius |5j . 
Second, in late stage spinodal decomposition, sponge- like 
network structures appear transiently in the asymmetric 
case 0, 0- The physical origin of these effects is now 
ascribed to the stress-diffusion coupling in viscoelastic 
binary mixtures This coupling is predicted to give 
rise to various effects including non-exponential decay in 
dynamic light scattering and shear-induced fluctua- 
tion enhancement. In the case of spinodal decomposition, 
the consequences of the dynamic coupling were studied 
analytically in the linear growth re gime |8« and numeri- 
cally in the late stage coarsening j^, Ei uM ■ 

Flow effects on phase separation in polymeric systems 
are even more dramatic Q, El E3- In this paper we 
consider a simple shear flow with mean velocity profile 



{v) = jye x , 



(1.1) 



where the flow is in the x direction, e x being the unit 
vector along the x axis, and the mean velocity gradient 
7 is in the y direction. Application of shear or exten- 
sional flow to viscoelastic systems in one-phase states 



sometimes induces a strong increase of the turbidity, 
indicating shear-induced composition heterogeneities or 
demixing. This is in marked contrast to shear-induced 
homogenization or mixing observed in near-critical flu- 
ids 0, and in ternary polymer mixtures 0, Il7| . 
In systems exhibiting shear-induced mixing, the entan- 
glement effects are not severe and the hydrodynamic in- 
teraction is suppressed by shear 0, Ei- In semidilute 
polymer solutions near the coexistence curve with high 
molecular weights (M >2 x 10 6 ), recent scattering exper- 
iments have most unambiguously detected shear-induced 
demixing d E3, H3, IH [23, HE H3, . Rheological ef- 
fects in sheared polymer solutions are also conspicuous, 
which include large stress fluctuations upon demixing by 
shear [2^, and a second overshoot in the shear stress 
as a function of time after application of shear [53, [2^ . 
Here microscope pictures of composition heterogeneities 
in polymer solutions and asymmetric polymer blends un- 
der shear are informative |29ll30j |. 

Theoretically, for sheared polymer solutions, relevance 
of the dynamical coupling was first pointed out by 
Helfand and Fredrickson [3l|. Some linear calculations 
for small fluctuations were also performed in Ginzburg- 
Landau schemes, where a conformation tensor represents 
the chain deformations [H |H 02 HH • Numerical anal- 
ysis (still in two dimensions) using such schemes gave in- 
sights into the nonlinear shear effects [Tiil l3rj l3lll3Sl l3^| . 
That is, slightly above the coexistence curve, composition 
heterogeneities 8<j> on mesoscopic spatial scales emerge 
with increasing 7. The amplitude of S(f> can even be of 
the order of the average (</>), while there are no clear in- 
terfaces. If this takes place, the system becomes turbid 
resulting in shear-induced phase separation or demixing 
observed above the coexistence curve. We remark that 
similar Ginzburg-Landau models 0> EH have been used 
to analyze shear-banding effects in wormlike micellar sys- 
tems mmi!^ 
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In entangled polymers the rheological relaxation time 
r can be very long |45|. so experiments in the New- 
tonian and non-Newtonian regimes are both possible, 
where the Deborach number De — 7T is larger or 
smaller than 1, respectively. In scmidilutc solutions, the 
strong composition-dependence of the solution viscosity 
rj = T]((f>) can lead to shear-induced fluctuation enhance- 
ment for weak shear 7T < 1 |3 1| , while the normal stress 
effect comes into play in the non-Newtonian regime 7T >1 

Without viscoelasticity a number of groups have per- 
formed simulations of phase separation in sheared simple 
fluids 0,E3- However, simulations have been still rare 
for sheared viscoelastic fluids [lj, |3g, |37j, |3a In 
this paper, we will examine nonlinear dynamic regimes 
of sheared polymer solutions in theta solvent below the 
coexistence curve. Both in simple and polymeric fluids, 
if the system is quenched below the coexistence curve un- 
der shear, the domain growth is eventually stopped and 
dynamically steady states are realized. For Newtonian 
fluids dynamics in sheared two-phase states have long 
been studied 0, H3, |H, , but for viscoelastic fluids 
such nonequilibrium effects remain almost unexplored. 

The organization of this paper is as follows. In Sec. 
II, we will explain our theoretical scheme. In Sec. Ill, 
we will present our numerical results in two dimensions 
below the coexistence curve for various shear rates and 
polymer volume fractions. Section IV summarizes new 
results and gives some predictions. 



II. THEORETICAL BACKGROUND 

In this section we briefly survey our theoretical frame- 
work to discuss viscoelastic phase separation under the 
shear flow (|f .1(1 . The gross variables in our model are 
the polymer volume fraction (f> and a conformation ten- 
sor W ■ The latter represents chain deformations and re- 
laxes slowly on the rheological time giving rise to a large 
viscoelastic stress. In Sec. II A we first present the vis- 
coelastic Gintzburg-Landau free energy, and then in Sec. 
II B the dynamic model is given. The rheological prop- 
erties of our dynamic equations are also discussed in Sec. 
II C. See Refs. 0, El for the details of our theoretical 
framework. 



where vq is the volume of a monomer (= a 3 with a being 
the monomer size in three dimensions) and \ is the inter- 
action parameter dependent on the temperature T (being 
equal to 1/2 at the theta condition). At the critical point 
we have = <f) c = A" 1 / 2 and A 1 / 2 (l - 2%) = -2. In the 
following it is convenient to scale </> and 2\ — 1 as 



A 1 / 2 (2 X -1). 



(2.2) 



In Fig.l we show the phase diagram in the plane of < 
and A 1 / 2 (l — 2%) = — u. The spinodal curve in Fig.l 
is obtained from (3 2 /fh/<90 2 )t = and is written as 
u = $ + 

To describe the viscoelastic effects on the composition 
inhomogeneities, it is convenient to introduce a tensor 

dynamic variable W— {W 7 ^ }, which is a symmetric ten- 
sor represen ting chain conformations undergoing defor- 
mations [3^, l.i.'St liiil IHol l5ll | . As shown in the following, 



the deviation 5Wi- 



Wi, 



Sij gives rises to a network 



stress. 

The Ginzburg-Landau free energy functional due to 
the fluctuations of <j> arL d W is written as 



F 



dr 



/fh + ^C|V0| 2 + - A GQ{W) 



(2.3) 



The coefficient C of the gradient term is calculated in 
the random phase approximation |4fij in the semidilute 
regime as 



C = {k B T/l%v )a 2 /(j>. 



(2.4) 



We shall see that G has the meaning of the shear modulus 
for small deformations changing rapidly compared with 
the rheological relaxation time r. It is assumed to be of 
the scaling form, 



G = {k B T/v ) g( j> o 



(2.5) 



where g is an important dimensionless parameter in our 
simulations, and is order of 1 in theta solvent [3 , Il3 | . 
Although experiments indicated a = 2.25 |52l l53j |. we 

will set a = 3 for simplicity. The simplest form for Q(w) 
is given by 



Q(w) - 5>w^ 



(2.6) 



A. Viscoelastic Ginzburg-Landau Free Energy 

Phase behavior of polymer solutions near the coexis- 
tence curve is usually described in term of the Flory- 
Huggins free energy density for the polymer volume frac- 
tion 6 assumed to be much smaller than 1 1451 . 



/fh = 



k B T 
vq 



(2.1) 



B. Dynamic Equations 

We next present the dynamic equations assuming the 
two-fluid dynamics for the polymer and the solvent with 
the velocities, v p and v s , respectively, including the new 

variable W [13, IH, l40| . These equations may be 
treated as Langevin equations with thermal noise terms, 
but we will neglect the noise terms hereafter. Formal 
frameworks for viscoelastic fluids have also been dis- 
cussed in the literature ISO, |5E • 
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We assume that the mass densities of pure polymer and 
solvent are the same. Then the polymer mass fraction 
Pp/p coincides with cf>. It obeys 
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-v • (<K) 

-V • (<bv) - V • \6(1 



w 



where 



(2.7) 



(2.8) 



is the relative velocity between polymer and solvent. The 
mean velocity v is defined by 



V = 0Vr 



(1- 



v s . 



(2.9) 



The two velocities v p and v s are expressed as v p = v + 
(1 — cf>)w and v s = v — <j)w. For simplicity, we assume the 
incompressibility condition for the average velocity, 



V • v = 0. 



(2.10) 



On the other hand, V • w is nonvanishing in general and 
gives rise to diffusion in (|2.7|) for small deviations around 
equilibrium. 

Because W represents the network deformation, its 
motion is determined by the polymer velocity v p and 
its simplest dynamic equation is of the form, 



d 



= — -SWi 



%3i 



(2.11) 



where Dij = 



dvpi/dxj is the gradient of the polymer 
velocity. The left hand side of (|2.11() is called the up- 
per convective time derivative in the rheological litera- 
ture [54[ , which keeps the frame invariance of the tensor 
properties of Wij . In l|2.21|l below, we will assume that 
the relaxation time r* on the right hand side depends on 
SWij as well as 4>. The usual rheological relaxation time 
r is obtained in the Newtonian limit: 



lim t* 



(2.12) 



In the problem of shear-banding flow, some authors re- 
placed 1/r* in (|?TT|) by (l/r*)(l - i 2 V 2 ) gflQ. 

The total stress tensor n= {IT-, } is expressed as 

(2.13) 



% = pSij +C(V i 0)(V i < 



where p is a pressure, Vi = d/dxi, a p = {<J P ij} is the 
network stress arising from the deviation SW^, and the 
last term is the viscous stress tensor with 770 being the 
solvent viscosity. Assuming low Reynolds number flows 



and setting dv/dt = (the Stokes approximation), we 
obtain [13| 



V- n= -Vpi + Tp + ri V 2 v = 0, 



(2.14) 



where pi = p — <pSF/S<p + /fh + C|V(/)| 2 /2 ensures the 
incompressibility condition (|2.1U|) . The T v is the force 
density acting on the polymer due to the fluctuation of 

<f> and W of the form, 



-<jN- 



.SF 1 



-QVG + V- <7p 



(2.15) 



The network stress tensor in l|2.13|l is expressed as 

a pij =G^2 W lk 8W kj + -GQ5 l3 . (2.16) 



To express v p and v s in terms of (j> and W, we as- 
sume slow processes and neglect the acceleration or in- 
ertial terms in the two- fluid dynamic equations 0. By 
setting dw/dt = 0, we obtain 



1 - 



(2.17) 



where £ is the friction coefficient between polymer and 
solvent and is estimated as £ ~ 67rryo^ 2 ~ rjQcf) 2 /a 2 in 
the semidilute solution in terms of the blob size The 
mean velocity v is expressed as 



v = (v) + 



1 



-r?oV^ p 



(2.18) 



where (v) is the mean flow such as the shear flow in 1)1.10 . 
[• • • ]j_ denotes taking the transverse part (whose Fourier 
component is perpendicular to the wave vector) and the 
inverse operation (— r/oV 2 ) -1 may be expressed in terms 
of the Oseen tensor in the limit of large system size • 
In the following, we make our eqations dimensionlcss 
by measuring space and time in units of I and tq defined 
by 



1 = 



aN 1 ' 2 



2V18 



1 _ 4fc B T 



iV 3 / 2 ' 



(2.19) 



where t is of the order of the gyration radius, and the 
time To is the conformation relaxation time of a single 
chain in the dilute case. In our simulations the mesh 
size in numerical integration will be set equal to I. The 
velocities will be measured in units of £/tq and the stress 
components given in (|2.13() will be measured in units of 



0q 



k B T/(v N 3 / 2 ) = tfc/(4r ). 



(2.20) 



To avoid cumbersome notation, in the follwoing, we use 
the same symbols for t, r, V, and the velocities even after 
rescaling. 
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C. Rheological quantities 

The conformation tensor W obeys H2.11|l with r* being 
replaced by t* /to in the dimensionless form. Following 
Refill we assume 

t*/t = (<S> 4 + 0.2)/(1 + Q), (2.21) 

where Q = Q(W) is defined by (|2.6I) and the factor 
1/(1 + Q) accounts for quickening of the stress relaxation 
under large deformations |55| . Similarly, some authors 
assumed a deformation-dependent stress relaxation time 
in the rheological constitutive equations 0, |5(| . 

Consequences of <|2.21[) are as follows, (i) In the dilute 
regime we have r* = 0.2tq. (ii) The relaxation time r in 
(12.12(1 in the linear response regime becomes 

t = t ($ 4 + 0.2). (2.22) 

(iii) The zero-frequency linear viscosity becomes 

77/770 = 1 + \g^r/T = 1 + iff* 7 , (2.23) 

where the small number 0.2 in r is omitted, (iv) Let us 
consider a homogeneous state under shear, where v p = 
v = jye x and (f> =const. In the high shear limit jr 3> 1, 
by solving i|2.11[l we obtain shear thinning behavior, 

°v*y ~ 5$ 3 (7r) 3/5 , 7V pl ~ ff <f 3 (7T) 4 / 5 . (2.24) 

Non-Newtonian behavior can arise from the factor 1/(1 + 
Q) in r* even in homogeneous states. 

Next we give dimensionless forms of the stress compo- 
nents for general inhomogeneous cases. From (|2.14|l and 
(|2.16|) the shear stress a xy and the normal stress differ- 
ence N\ = a xx — a vv in units of 00 are written as 

4 

dxy = (Jpxy ~ ^ V x $Vj,$ + 4(V a U„ + VyV x ), (2.25) 

+8(V x v x -V y Vy), (2.26) 
where a pxy and N p i are the network contributions, 

o pxy = g® 3 (W xx + W vy -l)W xy , (2.27) 
JVpi = g$ 3 (W xx + W vv -l)(W xx - W„„).(2.28) 

The second terms in <|2.25ll and (|2.26[) arise from inho- 
mogeneity in <I> and give rise to the surface tension con- 
tributions in two-phase states (5?]]. The last terms are 
the usual viscous contributions. 

The physical meaning of a pxy and AT p i can be seen if 
they are related to the degree of chain extension. To this 

end let us decompose SW as 

5W = wie 1 e 1 +w 2 e 2 e 2 , (2.29) 



where w\ and w 2 are the eigenvalues of SW with w\ >w 2 , 
and 

e\ = (cos 0, sin 0), e 2 = (— sin6*, cos 9) (2.30) 

are the corresponding eigenvectors with being the angle 
between the stretched direction and the x axis. We may 
assume — it/ 2 < 6 <tt/2 without loss of generality. In 
terms of these quantities we obtain 

°pxy = |$ 3 (l + wi +w 2 )(wi -w 2 ) sin 29, (2.31) 
N pl = ,g$ 3 (l + wi +w 2 )(w! - w 2 )cos29. (2.32) 

In weak shear rj <^ 1, we have 

wi — w 2 ~ T7, 9 — 7r/4 ~ r7, (2.33) 

so sin 29 = 1 and cos 29 ~ rj, leading to the well-known 
results a-o X V ~ 777 and N\ ~ 77T7 2 in the Newtonian 
regime |54j . Here, for the analysis in Sec. Ill, we in- 
troduce an extension vector defined by 

w(r,t) = (w 1 -w 2 )ei, (2.34) 

whose magnitude and direction represent the degree of 
chain extension and the extended direction. In our simu- 
lations we shall see that the magnitudes of w\ and w 2 are 
both considerably smaller than 1 at most space points. 
That is, the degree of extension is rather weak, but the 
network stress can overwhelm the viscous stress (cx 770) 
because of the large factor g& 3 . This justifies the Gaus- 
sian form of Q(W) in (|2.6|l in this work. 



III. NUMERICAL RESULTS 

We need numerical approach to understand the non- 
linear regime of shear-induced phase separation. To this 
end we integrate our model equations, (|2.7() and p.llfl . 
in the previous section on a 256 x 256 lattice in two di- 
mensions. The mesh size Ax is set equal to i in (|2.19|) 
[58| . We use a numerical scheme developed by one of the 
present authors 0> |H3 , which uses the deformed coordi- 
nates x' — x — jty and y' — y and enables the FFT (fast 
Fourier transform) method to be carried out for shear 
flow. Here we impose the periodic boundary condition 
f(x', y' + L) = f(x' + L, y') = f(x', y') for any quantity 
f(x',y') in terms of the deformed coordinates. 

As the initial condition at t = we assign Gaussian 
random numbers with variance (rms amplitude) 0.1 to $ 
at the lattice points. Then the initial $ at each point 
consists of the average (<&) and a random number. For 
(> we will solve the dynamic equations in the pres- 
ence of shear without the random noise terms. Here, if 
quenching is is not very far above the coexistence curve, 
the initial variance is not very small and our choice can 
be appropriate 60]. 
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A. Crossover from Newtonian to viscoelastic fluids 
below the coexistence curve 

Even below the coexistence curve, if the shear rate is 
sufficiently strong (yr ^> 1), the composition fluctuations 
vary in space gradually and there is no distinct phase sep- 
aration. However, with increasing the quench depth (or 
u) and/or decreasing the shear rate 7, the shear-induced 
composition fluctuations become composed of polymer- 
rich and solvent-rich regions with sharp interfaces. Here 
we show that the domain morphology strongly depends 
on the shear modulus G. 

In Fig. 2 we show phase separation patterns at t = 600 
and 2000 for g = (a), 0.01 (b), 0.1 (c), and 1 (d) after 
quenching at t = 0, where g represents the magnitude of 
G in (|2.5|l . The other parameters take common values 
as ($) = 2, u = 3, and 7 = 0.025. The initial Deborach 
number before phase separation is given by 7T = 7$ 4 = 
0.4. In the phase separated semidilute regions we have 

$cx = 3.70 (3.1) 

for u = 3 on the polymer-rich branch of the coexistence 
curve in Fig.l, so within the polymer-rich domains the 
viscosity ratio in l|2.23|l is given by 1 (a), 24.8 (b), 239 
(c), and 2384 (d). We can see gradual crossover from 
the patterns of the Newtonian fluids to those of highly 
asymmetric viscoelastic fluids. 

In Fig. 3 we show time evolution of the inverse of the 
perimeter length density for g =0, 0.01, 0.1, and 1, using 
the same conditions as in Fig. 2. It may be regarded as the 
typical domain size R(t). The domains continue to grow 
for g — and 0.01 within our simulation time t < 10 4 , 
but the growth is gradually slowed down with increasing 
g. For 3=1 and 0.1, we can see that the coarsening is 
nearly stopped. 



B. Domain size in two-phase flow 

Hereafter we fix g at 1. In Fig. 4 we show time evo- 
lution of the domain size R(t) for (<E>) = 2 and u = 3 
at various shear rates 7 =0.0005, 0.005, 0.025, and 0.05. 
Remarkably, the coarsening is faster for 7 =0.005 than 
for 7 =0.0005. This should be due to shear-induced 
coa gula tion of domains observed in near-critical fluids 
[T3L l6l|. where shear accelerates collision and fusion of 
the domains. For shear rates larger than 0.025, dynami- 
cal steady states are realized, where there should be a bal- 
ance between the thermodynamic instability and shear- 
induced domain breakup as in the Newtonian case [3.ll3|. 

The characteristic domain size i?D in the steady states 
are of interest. In Newtonian immiscible mixtures under 
shear we have Rp ~ C^ ew / (a xy ) in the low-Reynolds- 
number limit 0, Il3l llfjll. where CNew is of the order of 
the surface tension 7|62|- This formula follows from a 
balance between the surface energy density "f/Rr> and 
the shear stress. Also in our non-Newtonian case Fig. 5 



suggests the same form, 

i?D ~ G v { s /(a xy ) (3.2) 

at various 7. The coefficient C v i s is independent of 7 but 
dependent on (<J>) as 20, 13, and 11 for (4>) = 2.5, 2, and 
1.7, respectively. Fig. 5 displays the product (a xy )(t)R(t) 
as a function of time after quenching at g = 1 and u = 
3 with ($) = 2 and ($) = 2.5. The curves tend to 
composition-dependent constants independently of 7 at 
long times to confirm (|3.2|l . However, we cannot derive 
(|3.2|) using simple arguments in this paper. 

In the viscoelastic case, as in (|2.25|) and <|2.26[l . the 
stress consists of the three contributions. As shown in 
Fig. 6 at 7 = 0.05, the network stress a pxy dominates over 
the gradient contribution in the polymer-rich regions, 
whereas the network, gradient, and viscous ones are of 
the same order in the interface regions. The particularly 
large size of a pxy in the polymer-rich regions suggests 
that these regions should behave like percolated gels and 
mostly support the applied stress. On the other hand, 
the viscous contribution is very small in the polymer- 
rich regions, but is nearly the sole contribution in the 
solvent-rich regions. 

C. Fine domains at strong shear 

For the largest shear rate 7 = 0.05 in Fig. 4, Fig. 7 dis- 
plays time evolution of the domains, where the polymer- 
rich domains are percolated and the angle of extension 
6 defined by l|2~29"|> and (|2~3U|) is close to tt/4. These 
closely resemble the observed microscope pictures (in the 
xz plane) [HH3- The profile of $ at the bottom of Fig.7 
shows that $ becomes close to in the solvent-rich re- 
gions, while it is around $ cx in l|3.1|l in the polymer-rich 
regions. 

Fig. 8 displays the structure factor S(k x , k y ) of the com- 
position fluctuations in the steady state in Fig.7. It has 
sharp double peaks along the k x axis with peak wave 
number k p = (27r/256) x 9, although we cannot see 
marked anisotropy in the shapes of domains. The ori- 
gin of the peaks is that the domains are connected and 
hence are aligned perpendicularly to the flow direction 
on the spatial scale of 2tt jk v = 28. 

In addition, we notice that the profile of $ (bottom 
one in Fig.7) exhibit spike-like behavior around some ex- 
trema, where V$ varies over distances of order 5 [58| . 
These steep changes arise from convection due to the ve- 
locity fluctuations on such small scales in the nonlinear 
regime. 

Fig. 9 shows the average stress components (cr xy )(t) and 
(Ni)(t) vs time, which exhibit pronounced overshoots 
and subsequent noisy behavior. Here the network contri- 
butions in (|2.31l) and (|2.32|) are much larger than the sur- 
face tension contributions arising from the second terms 
in l|2.29|l and l|2.30|l by at least one order of magnitude. 
This is consistent with Fig. 6. Also shown is the average 
variance \J (5<& 2 )(t), which slowly increases over a long 
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transient time of order 2000. This arises from desorption 
of solvent from the polymer-rich regions into the solvent- 
regions, as was observed by Tanaka and coworkers 0,0- 



In fact Wolf and Sezen [64[ presented this view to in- 
terpret their finding of a viscosity decrease which signals 
onset of phase separation in shear in semidilute solutions. 



D. Large fluctuations at weak shear 

For the smaller shear rate 7 = 0.005, Fig. 10 demon- 
strates that the domain growth continues up to the sys- 
tem size in the simulation time t — 10 4 . Furthermore, 
comparing the two snapshots at t = 300 and 500, we can 
see that the domains are rotated as a whole in the early 
stage. Fig. 11 shows that the chaotic temporal fluctua- 
tions of the stress are much more exaggerated than in 
Fig. 9. Unusually the normal stress frequently takes neg- 
ative values, while it is always positive at 7 = 0.05 as 
shown in Fig. 9 at any ($). 

In Fig. 12 we further examine the origin of the strong 
fluctuations in this case. It displays the snapshots of the 
extension vector defined by 12.34f> . ((T xy ), (Ni), and the 
following rotationally invariant shear gradient 63] , 



) 2 - (v x uj, - v y v x y 



1/2 



(3.3) 



Notice 5 = for pure rotation. At t = 500, the an- 
gle 8 exceeds 7r/4 in most of the spatial points of the 
polymer-rich regions, resulting in (cos 29) = —0.21 and 
(Ni) = -0.44. At t = 1500, the points with 9 < tt/4 con- 
stitute a majority in the polymer-rich regions, leading to 
(Ni) = 0.24. These snapshots and those of the stress 
components clearly demonstrate the presence of stress 
lines forming networks, which are supported by the per- 
colated polymer-rich regions and where the extension and 
the stress take large values. The typical values of w\ — W2 
on these lines are 0.16 at t = 500 and 0.11 at t = 1500. 
For this shear rate we can see that (N\) becomes negative 
when the stress lines collectively rotate and the angle 9 
exceeds 7r/4 on these lines, while for 7 = 0.05 the stress 
lines are broken at much faster rates and (Ni) remains 
positive. On the other hand, S in Fig. 12 takes large 
values in the solvent-rich regions, indicating that the 
polymer-rich domains are largely rotated rather than be- 
ing anisotropically deformed in shear. This tendency be- 
comes conspicuous with decreasing shear. See the black 
regions of S in Fig. 12, where <!> is small and slipping is 
taking place. 

Fig. 13 shows a bird view of the velocity gradient 
dv x /dy (right) and the corresponding snapshot of <E> at 
t = 2000, where we choose 7 = 0.025, an intermediate 
shear rate in Fig. 4. The other parameters are the same as 
in Figs.4-12. We recognize that the solvent-rich regions 
support this velocity gradient. This means that strongly 
deformed solvent-rich domains act as a lubricant serv- 
ing to diminish the macroscopic stress or viscosity. This 
lubricant effect can be effective even for a small volume 
fraction of the solvent-rich domains because of the strong 
contrast of the viscoelastic properties in the two phases. 



E. Domain morphology and rheology 
for various ($) 

Next we show that the domain morphology also 
strongly depends on the average polymer volume frac- 
tion ($). Fig. 14 displays snapshots of $(r, i) at 7 = 0.05, 
u = 3, and g = 1 for three compositions, ($) = 1.2, 1.7, 
2.0, and 2.5 (from left). The time is t = 1500, where 
steady states are almost reached in all these cases. The 
profiles of $ in the x direction at y = 128 are also shown 
in the upper part. We can see that 4> = $ cx = 3.7 within 
the polymer-rich regions. With increasing ($), the col- 
lision frequency among the domains increases and the 
domain size decreases. In fact, i?D determined from the 
perimeter length is 13.75, 10.83, and 6.67 for ($) =1.7, 
2.0, and 2.5, respectively. Even for ($) = 1.2 the 
polymer-rich domains collide frequently and the domains 
shapes largely deviate from sphericity. In Fig. 15 we plot 
time evolutions of the average shear stress (a xy ) for these 
values of (<&). In the dynamical steady state, the normal- 
ized shear viscosity increase = v/Vo — 1 is given 
by 2, 4, 6.5 and 15 for ($) = 1.2, 1.7, 2.0 and 2.5, re- 
spectively. It is remarkable that the stress overshoot is 
nonexistent at small ($) and gradually develops with in- 
creasing ($). 



Deformation of sponge-like domains at small 
shear 



Finally we examine how the sponge- like domain struc- 
ture observed by Tanaka and coworkers 0, is deformed 
by shear flow. On the left of Fig. 16 we show one example 
of such a domain structure without shear for ($) = 1.2, 
u = 3, and g = 1, where the polymer-rich regions are per- 
colated despite relatively small ($) 0,0]. The right part 
is the result under shear 7 = 0.005 with the other param- 
eters being common, where the applied strain jt is 1, 2, 5, 
and 10 at t = 200, 400, 1000, and 2000 for the given snap- 
shots, respectively. For this weak shear rate, the initial 
stage of phase separation is not much different from the 
case without shear but the coarsening is quickened as in 
Fig. 4 for ($) = 2. Here the polymer-rich domains are 
gradually thickened but remain highly extended in the 
flow direction, while the domains in Fig. 14(a) are torn 
into pieces because of larger shear 7 = 0.05. Also in this 
case the average shear stress (a xy )(t) and normal stress 
difference (Ni(t)) exhibit chaotic behavior and the prod- 
uct (a xy )(t)R(t) tends to a constant (~ 2.4). Again as 
in Fig. 11, the negativity of the normal stress difference 
is conspicuous in the early stage (not shown here). 
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IV. SUMMARY AND CONCLUDING 
REMARKS 

Though performed in two dimensions, we have nu- 
merically solved the two-fluid dynamic model of sheared 
semidilute polymer solutions with theta solvent, where 
the chain deformations are represented by the confor- 
mation tensor Wij. The free energy density depends on 
Wij as well as the polymer volume fraction <f> as in l|2.3[) . 
In our simulations, the initial value of the composition 
is random, but the thermal noise terms in the dynamic 
equations are neglected for t > 0. We summarize our 
main results. 

(i) With varying the parameter g representing the 
magnitude of the shear modus, we have examined 
spinodal decomposition as in Figs. 2 and 3, which show 
the crossover of the domain growth from Newtonian to 
viscoelastic fluids. The domain growth is nearly stopped 
for g >0.1 at relatively large shear within our simulation 
time. 

(ii) The domain growth in spinodal decomposition has 
been examined at g = 1 with varying 7 in Fig. 4. The 
domain size i?D in steady states becomes finer with 
increasing 7 and, as demonstrated in Fig. 5, it satisfies 
the relation (|3.2() . Its theoretical derivation is not given 
in this paper. In Fig. 6 the network, gradient, and viscous 
stress contributions are compared, which indicates that 
the polymer-rich regions mostly support the applied 
stress. In Fig. 9 the variance ^/ (<5$ 2 )(t) continues to 
increase even after the saturation of the stress and the 
domain size, which arises from slow desorption of solvent 
from the polymer-rich regions 0, Q • 

(iii) Fig. 7 displays the composition patterns at relatively 
large shear 7 = 0.05. Its structure factor in Fig. 8 
suggests that the domains are more correlated along 
the shear-gradient direction rather than along the flow 
direction. This is consistent with the experiments 

mm in 113. 

(iv) At 7 = 0.05, <& fluctuates in the region |$ - ($)| 
between and <& cx in Fig. 7, exhibiting the spike-like 
behavior. Another notable difference is that the variance 
saturates rapidly above the coexistence curve 0, |3|J 
but slowly at a larger value in Fig. 9. 

(v) At much smaller shear 7 = 0.005, the domain size 
becomes larger as in Fig. 10 and the stress fluctuations 
look much enlarged as in Fig.ll. as compared to the 
stress curves in Fig. 9 for 7 = 0.05. For this small shear, 
percolated stress lines are formed in the polymer-rich 
regions as in Fig. 12. There, we also notice that the 
non-rotational velocity gradient S defined by (|3.3|) 
becomes large in the solvent-rich regions. Fig. 13 shows 
that the velocity gradient dv x /dy takes large values in 
the solvent-rich regions, suggesting that they act as a 
lubricant. 



(vi) Fig. 14 illustrates how the domain morphology 
depends on the average polymer volume fraction ($). 
Fig. 15 shows that the stress-strain curves also strongly 
depends on ($). 

(vii) Fig. 16 shows time evolution of the sponge- like 
patterns without shear d Q and under weak shear. 

Furthermore, we make comments on experimental as- 
pects. 

(i) We have calculated the space averages of the stress 
components. In usual rheology experiments, however, 
the stress acting on the boundary surface is measured. 
In a dynamical steady state, the time averages of the 
space-averaged stress and the surface stress do coincide, 
but their time-dependent fluctuations can be different. 
We note that the time-dependent fluctuations in these 
quantities become significant when the fluctuating enti- 
ties are of very large sizes. Here it is worth noting that 
more than three decades ago Lodge [26] reported abnor- 
mal temporal fluctuations of the normal stress difference 
at a hole of 1 mm diameter from a polymer solution con- 
tained in a cone-plate apparatus. He ascribed its origin 
to growth of inhomogeneities or gel-like particles. 

(ii) Below the coexistence curve, we have shown sharp 
stress overshoots at relatively strong shear in Fig. 9. How- 
ever, the stress-strain curves at relatively small shear in 
Fig.ll and those corresponding to Fig. 16 (not shown) are 
more complex, where the frequent negativity of the nor- 
mal stress is conspicuous. In the snapshots in Figs. 7, 10, 
and 16 the domains are collectively rotated as a whole in 
the early stage. Thus the first normal stress might be- 
come negative (at least) after its first peak. 

(iii) Disappearance of the stress peaks at small average 
volume fractions in Fig. 15 should also be observed. 

(iv) Experiments are needed on establishment of dy- 
namical steady states with the characteristic heterogene- 
ity length given by (|3.2() . Strong deviations of the 
polymer-rich domains from sphericity are also character- 
istic. These results gained from the simulations could be 
compared with microscope observations 0, |2t| . 

In future, we should perform simulations to examine het- 
erogeneous fluctuations in various sheared complex fluids 
such as polymer blends or wormlike micellar systems. 
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FIG. 1: Coexistence curve (solid line) and spinodal curve (dashed line) for polymer solutions obtained from 12. H in the plane 
of — u = A rl ^ 2 (l — 2\) and $ = 4 > / < t ) c- The points (x) represent the initial states of our simulations. 




FIG. 2: Crossover of domain patterns of $(r, t) below the spinodal curve with increasing g from to 1 for u = 3 and ($} = 2 
in shear flow with 7 = 0.025. 
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FIG. 4: Time evolution of the domain size R(t) below the spinodal curve for g = 1, u = 3, and (<&) = 2. Here 7 =0.0005, 
0.005, 0.025, and 0.05. At small shear rates flow-induced coagulation accelerates the domain growth as demonstrated by the 
curve of 7 = 0.005. For 7 ^> 0.005, shear-induced domain breakup becomes dominant and dynamical steady states are realized 
at smaller domain sizes. 
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FIG. 5: Time evolution of (a xy )(t)R(t) below the spinodal curve at various shear rates for g = 1 and u — 3. Here ($} = 2 
(solid lines) or ($} = 2.5 (dotted lines). At long times the curves tend to fluctuate around 15 — 20 confirming (3.4). 




FIG. 6: Three contributions, a pxy (network), — 4(V a; ( I>)(Vy ( I>)/ ( f> (gradient), and 4:(X7 x v y + V y v x ) (viscous) to the local shear 
stress a xy (x,y) in (2.40) in the range < y < 128 at x = for u = 3, ($) = 2, and 7 = 0.05 . The network stress is 
overwhelming in the polymer-rich regions, while the viscous one is relatively large in the solvent-rich regions. The gradient 
contribution is appreciable only in the interface regions, where the three contributions are of the same order. 
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FIG. 7: Time evolution of $(r, t) at 7 = 0.05 for w = 3, (<&) = 2, and g = 1 below the spinodal curve. The bottom figure is the 
profile at t = 500 in the x direction at y — 128. The domain size evolution is given by the curve of the largest shear 7 = 0.05 
in Fig. 4. 




FIG. 8: Structure factor S(k) for the composition patterns in Fig. 7 in the dynamical steady state. 
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FIG. 9: Time evolution of the average shear stress {a xy )(t), the average normal stress difference (Ni)(t), and the average 
variance \J (<5$ 2 }(t) (dotted line) at j — 0.05 for u = 3, ($) = 2, and g = 1 below the spinodal curve. 
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FIG. 10: Time evolution of $(r, t) at 7 = 0.005 for w = 3, ($} = 2, and g = 1 below the spinodal curve. For this shear rate 
Fig.4 shows that the domain size increases up to the system size at t = 10 4 . 
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FIG. 11: Chaotic time evolution of the average shear stress {(T X y)(t) and the average normal stress difference (Ni)(t) at 
7 = 0.005 in the run which produced Fig. 10. For this weak shear the deviatoric stress components exhibit large fluctuations 
and (Ni)(t) frequently takes negative values. 



<®>=2 u=3 g=1 Y=0.005 




FIG. 12: Patterns of the extension, the shear stress, the normal stress difference, and the rotationally invariant shear gradient 
13.31 at t — 500 (top) and t = 1500 (bottom). The corresponding composition patterns are shown in Fig. 10. We can see stress 
lines with large values of the extension and the deviatoric stress components. The S is large in solvent-rich slipping regions. 




FIG. 13: $(r,t) (left) and dv x (r,t)/dy (right) at j = 0.025 for u = 3, (*) = 2, and g = 1 in a steady state. This velocity 
gradient takes large values where slipping is taking place in the percolated solvent-rich region. 
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FIG. 14: Snapshots of $(r, t) for ($) = 1.2, 1.7, 2.0, and 2.5 in steady states at 7 = 0.05. The profiles of <E> in the x direction 
at y = 128 are also shown. 
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FIG. 15: Time evolution of {cr xy )(t) at various ($}. The corresponding composition patters are shown in Fig.14. The overshoot 
disappears at small ($}. 
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FIG. 16: Sponge-like domain structures without shear (left) and under shear 7 = 0.005 (right) for ($) = 1.2, u — 3, and g = 1. 



